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Abstract. The usual thermodynamic limit for systems of classical self-gravitating point particles 
becomes well defined, as a dynamical problem, using a simple physical prescription for the cal- 
culation of the force, equivalent to the so-called "Jeans' swindle". The relation of the resulting 
intrinsically out of equilibrium problem, of particles evolving from prescribed uniform initial con- 
ditions in an infinite space, to the one studied in current cosmological models (in an expanding 
universe) is explained. We then describe results of a numerical study of the dynamical evolution 
of such a system, starting from a simple class of infinite "shuffled lattice" initial conditions. The 
clustering, which develops in time starting from scales around the grid scale, is qualitatively very 
similar to that seen in cosmological simulations, which begin from lattices with applied correlated 
displacements and incorporate an expanding spatial background. From very soon after the formation 
of the first non-linear structures, a spatio-temporal scaling relation describes well the evolution of 
the two-point correlations. At larger times the dynamics of these correlations converges to what is 
termed "self-similar" evolution in cosmology, in which the time dependence in the scaling relation 
is specified entirely by that of the linearized fluid theory. We show how this statistical mechanical 
"toy model" can be useful in addressing various questions about these systems which are relevant in 
cosmology. Some of these questions are closely analagous to those currently studied in the literature 
on long range interactions, notably the relation of the evolution of the particle system to that in the 
Vlasov limit and the nature of approximately quasi-stationary states. 

INTRODUCTION 

The evolution of systems of a very large number of classical point particles interacting 
solely by Newtonian gravity is a paradigmatic problem for the statistical physics of 
long range interacting systems. It is also a very relevant limit for real physical problems 
studied in astrophysics and cosmology, ranging from the formation of galaxies to the 
evolution of the largest scale structures in the universe. Indeed current theories of the 
universe postulate that most of the particle-like clustering matter in the universe is 
"cold" and "dark", i.e. very non-relativistic and interacting essentially only by gravity, 
which means that the purely Newtonian approximation is very good one over a very 
great range of temporal and spatial scales. In practice, as we discuss further below, the 
very large numerical simulations used to follow the evolution of the structures in the 
universe, starting from the initial conditions on the density perturbations inferred from 
measurements of the fluctuations in the cosmic microwave background, are completely 
Newtonian. While simulations in this context have developed impressively in size and 
sophistication over the last three decades, the results they provide remain essentially 
phenomenological in the sense that our analytical understanding of them is very limited. 



Our first aim in this contribution is to clarify, in a statistical physics language, pre- 
cisely what the relevant problem is which is currently studied in the context of the prob- 
lem of structure formation in cosmology. The essential point is that it corresponds simply 
to a specific infinite volume limit of the Newtonian problem. A well defined limit of the 
resulting equations of motion is given by the case that the universe does not expand. This 
case, we explain, corresponds to the most natural definition of the usual thermodynamic 
limit of the Newtonian problem. We note that this limit is distinct from the "mean-field" 
or "dilute" one sometimes considered 1 in the literature, and in which, under certain 
circumstances, one can define thermodynamic equilibria. In the limit we consider the 
system has no thermodynamic equilibrium, and the system is always intrinsically out of 
equilibrium. 

The main body of this contribution will then be given over to the study of this (usual) 
thermodynamic limit of gravity, for a particular set of initial conditions. We will describe 
and characterise the evolution of clustering as observed in a set of numerical simulations, 
summarising the results of a series of recent publications with our collaborators [1, 2, 3]. 
Further we will highlight the fact that this evolution is qualitatively very similar to that 
seen in cosmological simulations (in an expanding universe). It is characterised notably 
by the fact that it is "hierarchical", which means that the non-linear structures form at 
a continually increasing scale, propagating towards larger scales at a rate which is fixed 
by the spectrum of initial fluctuations in the initial conditions. The evolution is observed 
also to tend asymptotically to what cosmologists call a "self-similar" behaviour, which 
means that the correlations manifest a simple spatio-temporal scaling behaviour. While 
the temporal part of this evolution may be explained by a simple linear fluid model for 
the growth of perturbations, the form of the spatial part remains, as in the cosmological 
model, unpredicted. 

We will then outline a model we have developed which describes very well the 
evolution of our simulations at early times. It is a two phase model, linking a perturbative 
treatment of the evolution at early times to an approximation in which the force on 
particles is dominated essentially by their nearest neighbour. 

In the last part we will turn to the relevance of our analysis to open unresolved ques- 
tions in the cosmological problem. Indeed one of our motivations has been that we may 
be able to use our system, which is relatively simple compared to a cosmological simu- 
lation, in this way as a kind of "toy model". We focus here specifically on the "problem 
of discreteness" in cosmological simulations: these simulations simulate a number of 
particles which is many orders of magnitude less than that in the physical problem, and 
the question is how this limits the resolution of the corresponding results. Stated more 
formally, the desired physical limit is the appropriate Vlasov-Poisson system, and the 
problem is to understand the finite N effects 2 . We discuss what we learn about this prob- 
lem from the case of the shuffled lattice. 

In our conclusions we suggest some directions for future work. In particular we 
discuss briefly the interest of relating out results to now very popular phenomenological 



See contribution of P.H. Chavanis. 
2 We will explain below that in the infinite volume limit it is, more precisely, a question of the relation 
between the evolution of the particle system, with finite particle density, and that in the Vlasov limit. 



models for the evolution of cosmological simulations, so-called "halo models". The 
central objects in this description (halos) have, apparently, very similar properties to 
those observed to be formed through the violent relaxation of finite self-gravitating 
systems. We note that the nature of these states of the finite system is analogous to 
the quasi- stationary states observed and much studied in one dimensional toy models 3 . 



NEWTONIAN GRAVITY IN THE INFINITE VOLUME LIMIT 

Let us start from a system of N identical point particles interacting by purely Newtonian 
gravity, i.e., with equations of motion 

where dots denote derivatives with respect to time, m is the mass of the particles, G is 
Newton's constant, and r ( - is the position of the z'-th particle. If iV is finite the evolution 
is well defined 4 , whether the system is open (i.e. in an infinite space), or enclosed in a 
region (with some appropriate boundary conditions). 
Let us consider now the usual thermodynamic limit, i.e., 

N 

N — > oo, V — > oo, at no = — = constant (2) 

For example we can consider taking N points which we randomly distribute in the 
volume V, increasing the number N in strict proportionality with the volume V. A little 
consideration of Eq. (1) shows that there is a problem when we take the limit: the sum 
on the right hand side, giving the gravitational acceleration of any given particle, is not 
unambiguously defined by this limiting procedure. To see this more explicitly, let us 
suppose that we choose our volume to be a spherical volume V(R S , To) centred on some 
arbitrary fixed point, at vector position To, and of radius R s , and that we take the infinite 
volume limit by increasing this radius R s — > oo. The force per unit mass for the point r, 
is then 

F(r,-) = -Gm lim £ ^'"V (3) 
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3 See the contribution of S. Ruffo in this volume. 

4 This is true of course assuming that r, ^ Yj for all ; ^ j. Starting from initial conditions in which this is 
true — as we always shall here — this undefinedness of the problem associated with this singularity of the 
force is never relevant, in a finite time. Whether this singularity is "left naked" or removed by a smoothing 
at some scale is, however, of course relevant to the very long time behaviour of the system. Indeed the 
so-called "gravothermal catastrophe" in these systems arises from the divergence of the microcanonical 
entropy due to configurations in which pairs of particles approach one another arbitrarily close, allowing 
other particles in the system to explore an unlimited phase space (in momentum). We will not explore the 
question here of the behaviour of these systems in the limit that t — > °°. See further comments below when 
we discuss the use of a regularisation at small scales of the force in numerical simulations. 



Defining w,-(r) = J^ ; y,-5(r — ry), and <5n,-(r) = J^ y y,-5(r — r,-) — no [where 5(r) is the 
Dirac delta function], this can be rewritten as 
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Force due to mean density 



The first term inside the brackets, representing the force per unit mass on the particle 
due to the mean density hq on a particle inside the sphere, gives, using Gauss' theorem, 
for any point r ( - e V(/??, r^) 



F„(rO = -^(r,-r„). 



(5) 



This result is independent of R s , and thus is also the expression of this contribution to 
the force in the limit that R s ; — > oo. We see, however, that it depends explicitly on the 
choice of origin. We will analyse the second term in the square brackets in Eq. (4), due 
to the correction to uniformity, in the next subsection. We will verify that this latter 
contribution may converge to a well defined value, independently of how the limit is 
taken, and thus that it cannot remove this dependence of the full force on the chosen 
point. 

To make the problem of Newtonian gravity in the usual thermodynamic limit well 
defined we need a uniquely defined force acting on particles. Given the arbitrariness in 
the result of the above calculation, we must thus give a prescription for calculating this 
force. As the ambiguity arises from the uniform component of the particle density, an 
evident choice for such a prescription is one which makes the contribution of the mean 
density zero for any particle. This is attained by any prescription for the sum in which 
the force on any given particle is determined by summing in a symmetric manner about 
that particle, e.g., 
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(7) 
(8) 



where V(R s ,Ti) is the spherical volume of radius R s centred on the i-th point. The force 
on a particle is then that due only to the fluctuations around the mean density. 

This prescription for the calculation of the force is essentially equivalent to what 
is known as the "Jeans' swindle" [4, 5]. To treat the dynamics of perturbations to an 
infinite self-gravitating pressureful fluid, Jeans "swindled" by perturbing about the state 



of constant mass density and zero velocity, which is not in fact a solution of the fluid 
equations he analysed. Formally the "swindle" involves removing the mean density from 
the Poisson equation, so that the gravitational potential is sourced only by fluctuations 
to this mean density. This is evidently equivalent to the prescription we have just given 
for the gravitational force, as seen explicitly in Eq. (8). 

A very nice discussion of the mathematical basis of this "Jeans' swindle" is given by 
Kiessling[6], who emphasises that the presentation of what Jeans did as a "mathematical 
swindle" is misleading. The "swindle" can in fact be given a perfectly firm mathematical 
basis when it is understood, as we have presented it above, as a physical regularisation 
of the force. The mathematical inconsistency in Jeans' analysis arises because it is done 
in terms of potentials which are badly defined, even in the thermodynamic limit as we 
have defined it. In Ref. [6] Kiessling shows that the Jeans swindle corresponds to a 
regularisation of the force given by the prescription used above in Eq. (6), which just 
corresponds to the physical prescription that the force on all particles be zero for an 
infinite constant density. Kiessling notes that an equivalent more physical form of the 
prescription is 

F(r f ) = -Gm lim V fi— % e -^ Ti -^ (9) 

where the sum now extends over the infinite space 5 , i.e., the value of the gravitational 
force in the thermodynamic limit is given as the limit of a screened gravitational force. 
One may thus actually consider the dynamics we describe below as characterising that 
of a system of particles with a screened gravitational force, but where the screening scale 
(= is much larger than the length scales on which "non-linear" clustering develops 
(see below). 



Force due to fluctuations 

Having given the prescription which makes the force zero in the usual thermodynamic 
limit for the case of an exactly uniform distribution of matter, we must consider also 
whether this prescription makes the force well defined in a distribution of particles, in 
which there are necessarily deviations from perfect uniformity, due to the discrete nature 
of the particles. To answer this question we must evidently consider how these particles 
are distributed in the infinite space, i.e., the nature of the fluctuations which they give rise 
to. For any distribution which one proposes to study as initial conditions of the infinite 
self-gravitating, one must ensure that the force is then well defined. 

To prescribe (and describe) how particles are distributed in infinite space we use the 
language of stochastic point processes, i.e., we consider the points to be distributed by 
a stochastic process, with certain properties. We assume only that these processes have 



This means that the thermodynamic limit has been taken before the limit in pi. 



only the following properties (which we require to employ them usefully here): 



• spatial ergodicity, i.e., the value of relevant quantities charactering their fluctua- 
tions/correlations in the infinite volume limit may be recovered by taking the en- 
semble average. 

• statistical homogeneity, i.e., averages are invariant under any rigid translation of 
the system. 

• uniformity, i.e., they have a well defined non-zero mean density, no = 
lim/v,y->oo(Af/V) = ("( x )) (where (...) denotes the ensemble average). 

In this framework the force on a particle, and also the force field (per unit mass) at 
any point in space, are stochastic variables. Let us consider the latter quantity 6 , denoting 
it f(r). It may be written as an integral over the stochastic variable n(r) as 



f(r) = —Gm lim 
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(10) 



The question of the well-definedness of the force may then be phrased as follows: is the 
probability distribution function (PDF) of the variable f, defined at any finite R s , well 
defined in the limit 7? s — > oo? An answer may be given as follows. Using the convolution 
theorem for Fourier transforms (FT) we can infer that 7 

|P(k)|ocjr 2 |5~n(k)| 2 (11) 

where the tilde denote the FT (and we have used the fact that the FT of x/|x| 3 is 
proportional to k/|k| 2 ). The square of the variance of any stochastic variable is, up to an 
appropriate normalisation in the volume 8 , its power spectrum (or structure factor). The 
integral of this quantity, it can be shown easily (see e.g. [7]) is proportional to the one 
point variance of the variable, from which is follows here that 

(f 2 ) oc I d 3 kk- z P(k) (12) 

where P(k) is the power spectrum for the density field n(r) of the point process. If this 
variance of the force is finite, the force PDF is necessarily well defined. This means 
that a sufficient condition for the force to be defined in the usual thermodynamic limit 
is just that fc _2 P(k) be integrable. Since we are interested here only in the possible 
divergences in all these quantities due to the long distance behaviour of the fluctuations 9 , 



We analyse here, for simplicity, the unconditioned force field, i.e., the force per unit mass at an arbitrary 
spatial point, rather than the force on an actual point of the distribution. These quantities will in general be 
different, notably for the small scale properties related to divergences of the force at zero separation. For 
the study of convergence/divergence of the force due to the long-range nature of the force, the quantities 
are equivalent. 

7 To do this FT, which is valid only in the sense of distributions, it is more appropriate to write the 
thermodynamic limit as in Eq. (9). 

8 The exact definition will be given in the next section. 

9 Note that we have already implicitly made this assumption in writing Eq. 10, as we have ignored the 
divergences when x falls on an occupied point if there is no regularisation of the force at zero separation. 



it is therefore sufficient that we require lim^o^(k) = 0. Thus it follows that the force 
is well defined in the infinite volume limit ifP(k) diverges more slowly at small k than 
k~ l . This can be compared with the condition of uniformity of the point process, which 
requires that P(k) be integrable, and thus only that it not diverge faster than k~ 3 as 
k — > 0. We can conclude that the force PDF is, at least, well defined for this sub-class of 
all uniform point processes. 

A Poissonian point process (i.e. particles thrown randomly in the infinite volume with 
any finite mean density no), which corresponds to a constant P(k), should thus give a 
well defined force in the thermodynamic limit. The PDF of the force on a particle for 
this case (with the same implicit regularisation we have described) was first calculated by 
Chandrasekhar in the forties [8], and gives the so-called Holtzmark distribution (see [7] 
for the full expression and a simplified derivation). The PDF is indeed well defined, and 
thus the dynamical problem in the thermodynamic limit with the given prescription. The 
variance is, in fact, infinite due to the small scale behaviour of the force: the probability 
at very large forces decays slowly (as a power-law) due to the unbounded growth of the 
interparticle force in configurations in which two particles are arbitrarily close. Imposing 
a cut-off regulating the divergence in the force at small separations, it is simple to 
determine (see [7]) how the variance of the force, which is then finite, diverges as this 
cut-off goes to zero. We will consider in detail below the case of "shuffled lattice" initial 
conditions. A detailed treatment of the statistics of the force field in these distributions 
is given in Ref. [9]. The power spectrum P(k) in this case (see below) is proportional to 
k 2 at small k, and it may be shown explicitly that the variance of the force is therefore 
finite in the defined thermodynamic limit, modulo divergences due to the small scale 
behaviour of the force. Indeed it is shown in Ref. [9] that the exact gravitational force 
PDF decays exponentially at large values of the field in distributions of this type in which 
there is zero probability of particle overlap, but as in the Poisson case (with infinite 
variance) if there is a non-zero probability for such configurations 10 . 



Relation to expanding universe case 

The principal conclusion of the preceding two sections is that there is a natural 
thermodynamic limit for the dynamical problem of self-gravitating particles in the 
Newtonian limit. The equations of motion of the particles are simply Eq. (1), where 
the sum is now prescribed to be taken symmetrically about each particle 11 . Provided 
the distribution of particles considered has fluctuations which decay sufficiently rapidly 
at large distances, this prescription, we have seen, gives a force on particles which is 
well defined 12 . Further, as we will see further below, the behaviour of fluctuations at 



An analogous treatment of generic power law interactions can be found in Ref. [10]. 

1 1 In practice this sum will be performed in the numerical simulations described below using the Ewald 
sum method. In this technique (see e.g. Ref. [11]) the long distance part of the sum is converted to a sum 
in reciprocal space, from which the badly defined term at k = 0, due to the mean density, is removed. 

12 The power spectrum P(k) is the FT of the reduced two point density-density correlation function, and 
thus the constraints on the sufficiently slow divergence at small k of the former correspond to constraints 



arbitrarily large scales is a property of the distribution which is not modified by the 
gravitational evolution, so the force will then be defined at all times if it is defined in the 
initial distribution 13 . 

The simplest way to understand the relation of this dynamical problem to that treated 
in cosmology is by comparing these equations of motion to those which apply in 
numerical simulations of structure formation in the universe. These are written (see e.g. 
Ref. [12]) 

x ; + 2^x i = -^£^3 (13) 
a(t) a 3 f- x, x, •< 

where x, are the positions (in "comoving coordinates") of the particles, and a(t) is the 
expansion rate of the cosmological model. The sum on the right-hand side, defining 
the force per unit mass on the particles, is evaluated in precisely the way we have just 
discussed, i.e., symmetrically, so that the force on each point is zero when we neglect 
the density fluctuations associated with the points 14 . The thermodynamic limit we have 
described for the purely Newtonian problem simply corresponds to d{t) = 0, i.e. to the 
formally non-expanding limit of the equations. 

The equations of motion Eq. (13) for the cosmological case, where a{t) ^ and a(t) 
obeys the correct equation derived from the full general relativistic equations, can in fact 
be derived also in a purely Newtonian framework, but using a limit different to the one 
we have considered. To see this let us return to Eq. (1), but with the sum now evaluated 
as given in Eq. (3), i.e. choosing a specific centre and summing the force in spheres 
centred on it. Taking now only the contribution to the force due to the mean density, i.e., 
neglecting the effect of deviations from uniformity, Eq. (5) gives the equation of motion 
for any particle 

AnGn Q m 

r/ = g (r/-r ). (14) 

To derive this formula we only used the fact that hq is well defined at any given 
instant, but did not assume that it was independent of time. It follows that we may seek, 
consistent with our assumptions, a solution of these equations describing a simple global 
dilatation of the whole (spherical) system, i.e. of the form 

r i (t)-r = R(t)[r i (t )-r ] (15) 



on the rapidity of decay of the latter at large separations. Alternatively one can also relate P(k) to the 
variance of mass in a volume, the latter being given by the integral of the former weighted by the FT of 
the window function describing the volume (see e.g. Ref. [7]). For P(k — > 0) ~ k", with n < 1, this gives, 
for example, that the variance in particle number in spheres of radius decays as R 3 ~". We have therefore 
shown above that the gravitational force has a finite variance due to fluctuations at large distances if the 
variance of these fluctuations grows more slowly with R than R 4 . 

13 The linearized fluid theory (see below) describes the evolution fluctuations at large scales, giving a 
simple ^-independent evolution for the power spectrum at small k. The latter thus maintains its initial 
functional form. 

14 In practice, just as in the simulations we will describe below, the sum is normally implemented 
numerically using an Ewald summation method. 



where R(to) = 1. In this case R(t) must then satisfy the equation 
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where we have defined p (t) = mno(t), the mean mass density at time t, and po 
mno(to). This equation may be integrated directly to give 



SkG rpo k 



where K is a constant of the motion, proportional to the conserved "energy": 
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(16) 
P(t ) = 

(17) 
(18) 



where Mq is a mass density proportional to the (constant) mass inside the sphere defined 
by any particle and the origin of the coordinates. The dynamics of the dilatation is 
clearly equivalent to that of a single particle in a central gravitational field, with K = 
determining the critical escape velocity at the initial time. 

These equations are identical to the Friedmann equations, derived within the frame- 
work of general relativity (see e.g. [13]) and which describe an expanding (or con- 
tracting) universe, for the case that the expansion (or contraction) is sourced only by 
non-relativistic matter. In practice all current standard-type cosmological models have 
in common that the universe is described very well by these equations (i.e. in which 
the observed expansion of the universe is driven only by non-relativistic matter 15 ) in the 
epoch at which structure formation develops. The constant R(t) and k acquire a different 
physical meaning to that ascribed to them in the Newtonian derivation — the former de- 
scribes stretching in time of the spatial sections of the four dimensional space-time and 
the latter the curvature of these spatial sections — but this difference in interpretation is 
of no relevance in the treatment of the problem of structure formation. 

To recover the full equations used in this context, i.e. Eq. (13), it suffices now to 
change coordinates, defining the comoving coordinates of any particle by 



Ti(t)-r = R(t)xi(t) 



Substituted in Eq. (1) this gives 
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(20) 

(21) 
(22) 



15 We will remark below on the addition of a "cosmological constant", which is now believed to cause an 
acceleration of the observed expansion at recent epochs. 



where here the sum over j e V(R s ,x) is over all the particles in the comoving sphere 
of radius R s (i.e. over the physical sphere of radius R(t)R s at time t) centred at the 
comoving coordinate position x. Making the identification of R(t) with a{t), this last 
form is precisely Eq. (13), with the prescription for the summation of the force made 
explicit. Note that we have used Eq. (16) in the passage between Eq. (20) and Eq. (21), 
and in passing from Eq. (21) to Eq. (22) we have assumed that the force due to the 
fluctuations (i.e. once the term due to the mean density is subtracted) is well defined 
independently of how the sum is taken. The criterion for this are identical to those 
discussed above. 

The only missing element in this derivation with respect to the full cosmological sim- 
ulations currently studied is the co-called "cosmological constant". This enters only as a 
modification of the Friedmann equation and thus of the function a{t) in the Eq. (13). We 
note that it in fact can be derived in a completely Newtonian framework, by generalising 
the above treatment with an additional classical repulsive force acting between all pairs 
of particles and linearly proportional to their vector separation. Because of its linearity it 
plays no role in the evolution of fluctuations, other than that encoded in the modification 
it causes of the uniform expansion. 

In summary the equations considered in cosmology to model formation of structure 
in the universe can be understood in a purely Newtonian framework as those describing 
the dynamics of particles, enclosed in a very large sphere and initially almost uniformly 
distributed with velocities close to those described by a homogeneous dilatation of this 
sphere (modelling the Hubble expansion). Formally we can consider the limit of the 
corresponding equations of motion in which there is no expansion. Physically this limit 
corresponds not to a very large static sphere — there is no static solution for such 
a configuration — but to a different regularisation of the gravitational force, which 
we have explained is naturally considered as the usual thermodynamic limit of the 
Newtonian problem of particles distributed uniformly in an infinite space. 

In this precise sense one can therefore study the problem of cosmological structure 
formation in the non-expanding limit. It is a natural starting point for a statistical 
mechanics approach to this problem, and as we will see one which already has most 
of the qualitative features of the cosmological problem. 

Relation to other N — > °° limits 

It is useful to comment on the relation of the (usual) thermodynamic limit we have dis- 
cussed to other limiting procedures used for the description of self-gravitating systems 
studied in the literature. More specifically we have in mind two such limits: 

• the "dilute" thermodynamic limit, and 

• the Vlasov-Poisson limit. 

The first of these limits arises in the context of the well-known treatment of gravity in 
the framework of the microcanonical ensemble, which has been discussed by a series of 
authors, going back to Emden in the early twentieth century (see e.g. [14] for a review, 
or the contribution of Chavanis in this volume). Such a treatment of a self-gravitating 



system of particles becomes possible if the particles are enclosed in a finite box and a 
cut-off is introduced to regularise the gravitational potential at the origin. Equilibrium 
states, i.e. states maximising the entropy exist for sufficiently large energy, which may 
then survive as meta-stable equilibria for the case that the cut-off goes to zero. This 
treatment becomes formally valid in the limit 

N 

N — > oo V — > °°, at — r-pr = constant. (23) 
yi/3 

The mean mass and number density thus scale as 

Po - n = ^ oc y- 2 /3 . (24) 

This scaling is clearly different to the limit we have taken, at fixed mean mass and 
number density. This means that the results obtained from this microcanonical treatment 
are not expected to be relevant to the dynamics we will study. This can be seen by noting 
that the characteristic time scale for the dynamical evolution of the systems we study 
(and which we will in fact use as our time unit below), is the dynamical time, defined as 

%n = 7=== ■ (25) 

This physical time scale 16 thus goes to infinity in the "dilute" thermodynamic limit, i.e. 
the dynamics we describe never takes place. Alternatively we can say that the effects 
described by this "dilute" limit are strictly relevant to the systems we study only on 
times scales which are infinitely long compared to those which we consider. Indeed the 
physics described by this microcanonical treatment is understood to be relaxational, 
characterised by time scales which diverge with N when expressed in terms of the 
dynamical time-scales. As we are studying already an infinite N limit, this relaxational 
time scale diverges, and any result from this treatment of a finite, closed, system cannot 
describe the physics of our infinite system, of which any finite subsystem is always 

17 

open . 

The Vlasov-Poisson limit, on the other hand, is a limit which has been studied for 
finite systems (see e.g. [15]). Its validity as an approximation to the full dynamics of an 
Af-body system has been rigorously demonstrated [16], for gravity with a regularisation 
of the singularity in the potential, when the limit Af — > °° is taken at fixed volume and 
a fixed time, with the scaling of the particle mass m« \/N. This latter scaling makes 
the acceleration given by the sum on the right hand side in Eq. (1) converge to the 
appropriate mean field value. In this case we therefore have 

N mN 
uq = — and po = = constant . (26) 



16 E.g. for a mass density of one gram per cm 3 , one finds Tdyn ~ 18.2 minutes. 

17 This treatment may, nevertheless, be relevant, qualitatively, to understanding aspects of the very long 
time evolution of the systems we study. Indeed, as we will describe, highly clustered regions can behave 
over much longer timescales than Td yn as roughly independent from the rest of the system, and so may 
manifest long-time relaxational behaviors like those of isolated finite systems which are, at least partially, 
described by this microcanonical treatment. 



Again this is clearly a different limit to the one we have considered. We see, however, 
that this limit leaves invariant the characteristic time scale Td yn , which we noted is that 
of the dynamics for the infinite systems we will describe and analyse below. The reason 
for this is in fact that this time scale emerges from a treatment of this dynamics which 
is derived by treating it in a Vlasov-Poisson approximation. In this case, however, this 
approximation now applies to the infinite system, i.e., in which the limit N — > °° has 
already been taken, but as we have described, at constant no and po- Although there 
is no rigorous demonstration in the literature of the validity of such a limit for such a 
system, it is clear that it must correspond, if it indeed exists, to the limit 

n — > oo at po = mno = constant . (27) 

Indeed when the Vlasov limit is taken in a fixed (finite) volume it may be written in 
this volume-independent way, which thus generalizes without difficulty to the infinite 
volume case. As we will discuss further below, in cosmology this limit is of great im- 
portance, as the aim of the numerical simulations of self-gravitating Newtonian particle 
systems in this context is to reproduce this limit as well as possible. 

Lastly we emphasize an important point concerning the usual thermodynamic limit 
for self-gravitating particles: we have shown that this limit exists for the dynamical 
problem. It is well known, and indeed one of the starting points of any discussion of 
long-range interactions, that in this limit the use of the usual instruments of equilibrium 
thermodynamics is problematic. Indeed the gravitational potential violates both the es- 
sential conditions necessary for their use formulated by Ruelle[17] — the condition of 
"temperedness" due to the slow decay of the potential (i.e. slower than the dimension- 
ality of space) at large separations, and the condition of stability due to the behaviour 
at small separations. Our treatment of the problem in the usual thermodynamic limit is 
completely consistent with this. Indeed it turns that the dynamical problem we have for- 
mulated is an intrinsically out of equilibrium one, the evolution of the system remaining 
always explicitly time dependent (in a calculable manner, determined as we will see, 
by the analysis in the Vlasov-Poisson limit). Thus any treatment of the problem as an 
equilibrium one, or a quasi-equilibrium one (see e.g. the contribution of Saslaw in this 
volume) will at best be a useful, but non-rigorous, approximation. 

EVOLUTION FROM SHUFFLED LATTICE INITIAL 
CONDITIONS: PHENOMENOLOGY 

We describe in this section, in a phenomenological manner, the evolution under their 
self-gravity of particles initially placed in a shuffled lattice (SL) configuration as ob- 
served in a set of numerical simulations. This synthesizes the salient results of a more 
detailed study reported in [1]. Having given the exact definition of a SL, and explained 
our choice of this specific class of initial conditions, we will present the results of sim- 
ulations, focusing on the standard two point characterisations of the observed evolving 
correlations, in terms of the reduced two point correlation function and its reciprocal 



space equivalent, the power spectrum . In the presentation of these results we will 
highlight the qualitative similarity with the known behaviour of analogous cosmolog- 
ical simulations, and explain for uninitiated readers the essential theoretical results — 
just the original analysis of Jeans of the evolution of density perturbations in the fluid 
limit — used to understand them. 



Shuffled Lattices 

We use the term SL to refer to the infinite point distribution obtained by randomly 
perturbing a perfect cubic lattice: each particle on this lattice, of lattice spacing £, is 
moved randomly ("shuffled") about its lattice site independently of all the others. A 
particle initially at the lattice site R is then at x(R) = R + u(R), where the random 
vectors u(R) are specified by p(u), the PDF for the displacement of a single particle. 

An ensemble of SL configurations is thus characterised by the lattice spacing £ 
and the PDF p(u). We will use in the simulations below a simple top-hat form, i.e., 
p(u) = (2A)~ 3 for u G [—A, A] 3 , and p(u) = otherwise 19 . Taking A -> 0, at fixed £, 
one thus obtains a perfect lattice, while taking A — > <*> at fixed £, gives an uncorrected 
Poisson particle distribution^]. Once the form of the PDF is specified the class of SL 
we consider is thus a two parameter family, characterized by £ and A. We refer to the 
latter as the shuffling parameter. Using the top-hat PDF it is simple to verify that it is 
simply the square root of the variance of the displacement PDF, i.e., A 2 = / d 3 uu 2 p(u). 
It is convenient to define also the dimensionless ratio 6 = A/£, which we will refer to as 
the normalized shuffling parameter. 

To characterize the correlation properties of the SL let us consider its power spectrum. 
We recall that for a point (or continuous mass) distribution in a cube of side L, with 
periodic boundary conditions, this quantity is defined as 20 

P(k) = l(|5(k)| 2 ), (28) 

where 8 (Is) are the Fourier components of the density contrast 8(x) = (p(x) — Po)/po 
(and p(x) is the mass density), i.e., 

8(k)= I 5(x)exp(-/k-x) d 3 x. (29) 
Jl 3 

for k = (2n/L)n where n is a vector of integers. The infinite volume limit is then 
obtained by sending L — > oo 5 keeping the average mass (or number) density fixed. We 
give the finite volume form for P(k) here as we will in fact treat in our simulations an 
infinite system defined in this way. 



We will recall the definitions of these quantities and some of their essential properties at the appropriate 
points, but not the details of how they are estimated in the simulations. These may be found in [1]. 

19 We will explain below that the exact form of the PDF is not of importance. 

20 We adopt the normalisation which is used canonically in cosmology, for which the asymptotic be- 
haviour at large k, characteristic of any point process, is P(k — > °°) = 



It is straightforward (see Ref. [7]) to show that the power spectrum for a generic SL 
is given by 

P(k) = 1 -I^ k )| 2 +L 3 £ |# k )|2fc( M ^) (30) 
no n ^ I 

where p(k) is the Fourier transform of the PDF for the displacements p(u) (i.e. its 
characteristic function), and 5k is the three-dimensional Kronecker symbol. The second 
term in Eq. (30) is non-zero only at the corresponding discrete (non-zero) values of k. 
It is in fact simply the power spectrum of the unperturbed lattice modulated by |p(k) | 2 . 
Thus only the first term in Eq. (30) contributes around k = 0. Taylor expanding this term, 
assuming only that p(u) has finite variance, we obtain the leading small-fc behavior 

P( k ) « k ^L = IfcW, (31) 

where k = |k| . Note that this small-fc behavior of the power spectrum of the SL therefore 
does not depend on the details of the chosen PDF for the displacements, but only on its 
(finite) variance. For this reason in fact the results we will discuss for the chosen PDF 
do not in fact depend in detail on this specific form 21 . Note that the power spectrum is a 
function only of k at small k, which means that the large scale fluctuations in the system 
are statistically isotropic. The behavior P(k) k 2 means that, as one would expect, an 
SL is at large scales a distribution with fluctuations which are suppressed compared to 
those in a Poisson distribution. Indeed the SL belongs for this reason to the class of 
"superhomogeneous"[18] (or "hyperuniform"[19]) point processes with this property, 
characterized by the behavior P(k — > 0) = 0. 

Let us now comment on why we choose this particular class of SL initial conditions. 
The reason is that they are a very simple (two parameter) class of initial conditions 
resembling those used in cosmological simulations. In this context initial conditions are 
prepared by applying correlated displacements to a lattice. By doing so one can produce, 
to a good approximation, a desired power spectrum at small wavenumbers 22 . The spectra 
of currently favored models contain several parameters, roughly interpolating between 
power law forms P(k) k n , with n varying from +1 at the smallest k to close to —3 at 
the largest k included in the initial conditions of numerical simulations. In particular the 
SL allows us to mimic an important feature of these simulations, related to the problem 
of discreteness we will discuss further below: with the two parameters characterising 
the class of initial conditions, we can modify independently the particle density and 
the amplitude of the power spectrum at small wavenumbers. This can be seen from 
Eq. (31), as it is clear that to keep the latter fixed it suffices to impose, when varying £, 
the condition 

8 2 £ 5 = constant. (32) 



In Ref. [2] simulations like those discussed here, but with a Gaussian PDF are analysed. 
See [20] for a detailed analysis of the algorithm used. 



TABLE 1. Details of the SL used as initial conditions in the sim- 
ulations of [1]. N is the number of particles, L is the box size, I the 
lattice constant and A (8) the (normalized) shuffling parameter. 



Name 




L 


I 


A 


8 


m/ntfA 


SL64 


64 


1 


0.015625 


0.015625 


1 


1 


SL32 


32 


1 


0.03125 


0.0553 


0.177 


8 


SL24 


24 


1 


0.041667 


0.00359 


0.0861 


18.96 


SL16 


16 


1 


0.0625 


0.00195 


0.03125 


64 


SL128 


128 


2 


0.015625 


0.015625 


1 


1 



Numerical Simulations 

In numerical simulations we cannot treat of course a truly infinite SL. Instead we treat 
the infinite system defined by such a lattice, in a finite cubic box of size L = N^ 3 £, plus 
periodic copies. Our results will then be representative of the thermodynamic limit only 
insofar as they do not depend on the size of the box. We will see below that the nature 
of the dynamics is such that this can be true only for a finite time, i.e., the evolution of 
this infinite, but periodic, system will represent that in the thermodynamic limit for a 
finite time. To follow the evolution for a longer time requires an ever larger box. Any 
effects which depend on the box size are finite size effects which are not related to the 
real behavior of the well defined physical limit we are studying 23 . 

In Table 1 are given the details of the SL initial conditions of the N body simulations 
reported in [1], and which we will discuss here. 

The numerical evolution has been performed using the publicly available code GAD- 
GET [21, 22]. This is based on a tree algorithm for the calculation of the force, and 
allows one to perform simulations in an infinite periodic distribution using the Ewald 
summation method. The potential used is exactly equal to the Newtonian potential for 
separations greater than a softening length £, and regularized below this scale 24 . More 
details of tests on the accuracy of the simulations (energy conservation, robustness to 
change of integration parameters etc.) can be found in [1]. 

Let us explain the rationale for the choice of the parameters in Table 1. Firstly, we 
have chosen our (arbitrary) units of length, mass and time as follows. Our unit of length 
is given by the box side of the SL64 simulation and our unit of mass by the particle 
mass in this same simulation. The particle masses are then chosen so that the mass 
density po is constant, which means that the dynamical time defined above in Eq. (25) 
is the same in all simulations. This is convenient because this is, as we have mentioned 
above, an appropriate characteristic time-scale for the evolution which can be derived 
from the fluid limit (see below). We have chosen a smoothing parameter e = 0.00175 
in all simulations. It is therefore in all cases significantly smaller — at least a factor of 



J A good analogy is in the simulation of the out of equilibrium dynamics of glassy systems, e.g., the 
ordering dynamics of a quenched ferromagnet. 

We do not give here the rather complicated functional form of this smoothing, which can be found 
in [21, 22]. The important point here is that the smoothed force goes to zero continuously as particles 
approach one another. 
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FIGURE 1. The power spectra (averaged in spherical shells) of the SL configurations specified above 
in Tab. 1 as a function of the modulus of k. The solid line is the theoretical (°< k 2 ) behavior for small k 
given by Eq. (31). At large k, the four power spectra are equal to 1/riQ, with the corresponding value of 
no- The peaks arise from the second term in Eq. (30). From Ref. [1]. 



ten — than the lattice spacing (which is close in all cases to the initial average distance 
between nearest neighbors). The box size is the same in all but one simulation. This latter 
simulation (SL128), which is the biggest one, is used to test the accuracy with which our 
results are representative of the thermodynamic limit. Thus it is chosen to have the same 
parameters to SL64, differing only in its volume (which increases by a factor of 8). 
Finally the values of 8 have been chosen for each £ so that the condition in Eq. (32) is 
fulfilled. The measured power spectra of realizations of SL with the parameters given 
in Table 1 are shown in Fig. 1. We see, up to statistical fluctuations, that the spectra 
are indeed of the same amplitude at small k. Note that the curve corresponding to each 
initial condition is most easily identified by the fact that the asymptotic (large k) level 
of the power spectrum is inversely proportional to the mean particle density. 

Lastly, the particles are assigned zero velocity in the initial conditions (at t — 0), and 
they are run in each case until a time when the finite size of the box manifestly becomes 
important (see below). 

Let us comment further on the use of the smoothing parameter e. The goal here is to 
study the evolution of a distribution of self-gravitating point particles, and the reason for 
the introduction of £ is that it allows a much more rapid numerical integration. It does 
this because it modifies the trajectories in which particles have close encounters, which 
in lead to very large accelerations (and thus the necessity for very small time steps). 
The use of the smoothing is only justified therefore if the associated physical effects 
are not important for the macroscopic properties we will study. A priori we do not know 
whether this is the case, and we simply test numerically (see [1]) for the robustness of our 
results to changes, towards significantly smaller values, of the smoothing parameters. We 
interpret the observed £-independence of our results as meaning that the artificial effects 
it introduces in the dynamics may be relevant to the quantities we measure only on time 
scales longer than those we can study in our simulations. 

We mention in this context the following point, to which we will return further 
below. In cosmological simulations of particle systems the smoothing is understood 
as a physical parameter, introduced with the intention of making the particle system 




FIGURE 2. Snapshots of SL64, at t = 0, and the evolved configurations obtained at subsequent times, 
t=3,6,8. These are projections on the z — plane. 



behave more like the desired Vlasov-Poisson limit (in which two body encounters 
are neglected). In practice, however, values of the ratio e/£ even smaller than those 
adapted here are typically used. Thus the parameter choices we have made here would 
be considered as perfectly reasonable if our goal were, as in cosmology, to reproduce 
the Vlasov-Poisson limit (and, as far as possible, nothing else). 

Results 

In Fig. 2 are shown four snapshots 25 of the simulation SL64. We see visually that non- 
linear structures (i.e. regions of strong clustering) appear to develop first at small scales, 
and then propagate progressively to larger scales. Eventually the size of the structures 
become comparable to the box-size. From this time on the evolution of the system will 
no longer be representative of the thermodynamic limit we are studying. Up to close to 
this time, however, it is indeed the case that all the properties we will study below show 
negligible dependence on the box size (see Refs. [1,2] for more detail). 



The results taken from Ref. [1] are given in units which are not exactly those of dynamical time, but 
rather in which Td yn = 1.092. This corresponds to time units of 1000 seconds if the mass density were 
lg.cm~ 3 . 



Evolution of correlations in real space 



We consider now the evolution of clustering in real space, as characterized by the 
reduced correlation function t, (r) . We recall that this is simply defined as 

§(r) = (5(x + r)5(x)>, (33) 

where (...) is an ensemble average, i.e., an average over all possible realizations of the 
system (and we have assumed statistical homogeneity). It is useful to note that this can 
be written, for r ^ 0, and averaging over spherical shells, 

{W= W>£-1, (34) 
no 

where (n(r)) p is the conditional average density, i.e., the (ensemble) average density of 
points in an infinitesimal shell at distance r from a point of the distribution. Thus t, (r) 
measures clustering by telling us how the density at a distance from a point is affected, 
in an average sense, by the fact that this point is occupied. In distributions which are 
statistically homogeneous the power spectrum P(k) and | (r) are a Fourier conjugate 
pair (see e.g. Ref. [7]). 

The correlation functions we will discuss below will invariably be monotonically 
decreasing functions of r. It is then natural to define the scale X by 

£(A) = 1. (35) 

This scale then separates the regime of weak correlations (i.e. ^(r)cl) from the regime 
of strong correlations (i.e. £ (r) 1). In the context of gravity these correspond, approx- 
imately, to what are referred to as the linear and non-linear regimes, as a linearized treat- 
ment of the evolution of density fluctuations (see below) is expected to be valid in the 
former case. Eq. (35) can also clearly be considered as a definition of the homogeneity 
scale of the system. Physically it gives then the typical size of strongly clustered regions. 

In Fig. 3 is shown the evolution of the absolute value |£ (r)\ in a log-log plot, for the 
SL128 simulation. These results translate quantitatively the visual impression gained 
above. More specifically we observe that: 

• Starting from ^(r)<l everywhere, non-linear correlations (i.e. £ (r) ^> 1 ) develop 
first at scales smaller than the initial inter-particle distance. 

• After two dynamical times the clustering develops little at scales below e. The clus- 
tering around and below this scale is characterized by an approximate "plateau". 
This corresponds to the resolution limit imposed by the chosen smoothing. 

• At scales larger than £ the correlations grow continuously in time at all scales, with 
the scale of non-linearity [which can be defined, as discussed above, by t, (A) = 1] 
moving to larger scales. 

From Fig. 3 it appears that, once significant non-linear correlations are formed, the 
evolution of the correlation function t, (r) can be described, approximately, by a simple 
"translation" in time. This suggests that £(r,/) may satisfy in this regime a spatio- 
temporal scaling relation: 

$(r,0«S(r/tf 5 (0), (36) 
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FIGURE 3. Behavior of the absolute value of the correlation function \E,(r)\ in SL128 at times t — 
0, 1,2,4,6,8. The vertical dotted line indicates e. From [1]. 




8 (w. resc.) ^ 



10" 5 10" 4 10" 3 10" 2 10"' 10° io' 

r 

FIGURE 4. Collapse plot of E, (r, t ) : for each time t > 1 we have rescaled the x-axis by a time-dependent 
factor to collapse all the curves (dashed ones) to that at time t = 1. We have added for comparison 
E, (r,t = 8) without rescaling ("w. resc", continuous line). From [1]. 



where R s (t) is a time dependent length scale which we discuss in what follows. In order 
to see how well such an ansatz describes the evolution, we show in Fig. 4 an appropriate 
"collapse plot": ^ (r,t) at different times is represented with a rescaling of the x-axis by 
a (time-dependent) factor chosen to superimpose it as closely as possible over itself at 
t = 1, which is the time from which the "translation" appears to first become a good 
approximation. We can conclude clearly from Fig. 4 that the relation Eq. (36) indeed 
describes very well the evolution, down to separations of order £, and up to scales at 
which the noise dominates the estimator (see [1] for further details). In Fig. 5 is shown 
the evolution of the rescaling factor R s (t) found in constructing Fig. 4, as a function 
of time [with the choice R s (l) = 1]. Also shown is a theoretical curve which we will 
explain in the next subsection. Note that since we have defined the homogeneity scale X 
by £ (A) = 1 it is clear that, once the spatio-temporal scaling relation is valid, we have 
X(t) oc R s (t). A fit to a simple functional form for S(r), a shallow power law followed 
by an exponential cut-off, may be found in Ref. [1]. 
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FIGURE 5. Evolution of the function R s (t) in SL128 (points) compared with its prediction from 
linearized fluid treatment, as explained in the next section. From [1]. 

Evolution of the Power Spectrum 

To see how the observed temporal dependence of the correlation function is explained 
theoretically, in analogy to how this is done in cosmological simulations, we outline very 
briefly the essential result of the Jeans-type analysis of the evolution of perturbations 
in a self-gravitating, and pressureless, fluid. For comparison of its predictions with 
the numerical results for the discrete N-body system, it is simplest then to turn to the 
reciprocal space description of the clustering 26 . 

The equations describing the evolution of a non-relativistic pressureless self- 
gravitating fluid are the following 27 

d f p + V x -(pv) = 0, d fV + (v-V x )v = g 
V x -g = -4xG(p-p ), V x xg = 0, 

where p(x,f) is the mass density, v(x,?) the velocity field and g(x,?) the gravitational 
field. These equations can also be obtained by an appropriate truncation of Vlasov- 
Poisson equations. Linearizing in the density perturbations, and the velocity, one obtains 

8 (x, t) = AkGPqS (x, t) and 8 (k, t) = 4nGp 8 (k, t) . (37) 



26 The real and reciprocal space description are, of course, fundamentally equivalent. The convenience of 
one space over another is simply in how effects due to their estimation in finite samples enter. 

27 As mentioned above a rigorous derivation of these equations in the thermodynamic limit as we have 
defined it can be found in Ref. [6]. 
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FIGURE 6. Evolution of the power spectrum in SL128 (solid lines — label FG): the curves are for time 
equal to 0,2,4,6,8 (from bottom to up). The dashed lines labeled with LT show the predictions of fluid 
linear theory, i.e., Eq. (39) with P(k,0) measured in the simulation at t = for the same time steps. The 
arrow labeled shows the value of the corresponding Nyquist frequency — n/£. From [1]. 

For the case we are analysing, of zero initial velocities, therefore 

5(M) = S(k,0)cosh (^/ATtGpQ /) . (38) 

and thus for the power spectrum 

P(k, t) = P(k, 0) cosh 2 (t/r dyn ) . (39) 

The evolution of the power spectrum as estimated in SL128 is shown in Fig. 6. Along 
with the numerical results is shown the prediction Eq. (39). We observe that: 

• The linear theory prediction describes the evolution very accurately in a range 
k < k*(t), where k*(t) is a wave-number which decreases as a function of time. 
This is precisely the qualitative behavior one would anticipate (and also observed 
in cosmological simulations): linear theory is expected to hold approximately only 
above a scale in real space, and therefore up to some corresponding wavenumber 
in reciprocal scale, at which the averaged density fluctuations are sufficiently small 
so that the linear approximation may be made 28 . This scale in real space, as we 
have seen, clearly increases with time, and thus in reciprocal space decreases with 
time. We note that at t = 6 only the very smallest fc-modes in the box are still in this 
linear regime, while at t = 8 this is no longer true (and therefore finite size effects 
are expected to begin to play an important role at this time). 

• At very large wave-numbers, above k^, the power spectrum remains equal to its 
initial value l/«o- This is simply a reflection of the necessary presence of shot 
noise fluctuations at small scales due to the particle nature of the distribution. 




A more precise study of the validity of the linearized approximation is given in [1]. 



• In the intermediate range of k the evolution is quite different, and slower, than that 
given by linear theory. This is the regime of non-linear clustering as it manifests 
itself in reciprocal space. 

Spatio-temporal scaling and "self -similarity" 

We now return to real space and explain a little further how the spatio-temporal scaling 
relation Eq. (36) for the correlation function we have observed may be linked to the 
linearized fluid theory, and then what it tells us about the nature of the clustering in the 
system. 

In the context of cosmological N body simulations this kind of behavior, when R s {t) 
is itself a power law (in time), is referred to as self -similarity. Such behavior is expected 
in an evolving self-gravitating fluid (see e.g. [13, 23, 24]) because of the scale-free 
nature of gravity, if the expanding universe model and the initial conditions contain 
no characteristic scales, i.e., if the input model has a simple power law P(k) °< k". On 
theoretical grounds there are different expectations ([24, 25, 23]) about the range of 
exponents n of the power spectrum which should give self-similar behavior, because 
of how effects coming notably from particle discreteness, the finite box size and force 
smoothing break the scale-free idealisation. It is widely agreed that it applies well for 
— \ < n < — lin the literature, but there are differing conclusions about the observed 
degree of self- similarity outside this range. 

The simple arguments used in this context to predict the temporal scaling of the 
correlation functions can be generalized easily to the case of a non-expanding universe 
(which, just like the relevant cosmological models, has no characteristic scale). To do so 
one simply assumes that such a spatio-temporal scaling relation holds exactly, i.e., at all 
scales, from, say, a time t s > 0. For t > t s we have then 

P(k,t) = [ exp(-/k-r)^(r,0d 3 r (40) 

Jl 3 

= R 3 s (t) [ exp(-^,(0k-x)S(|x|)d 3 x (41) 

= R 3 s (t)P(R s (t)k,t s ) . (42) 

where we have chosen R s (t s ) = 1. Assuming now that the power spectrum at small k is 
amplified as given by linear theory, i.e., as in Eq. (39), one infers for any power spectrum 
P(k) ~ k n that 29 

2 

m = V^r¥] "* - «p \ n 2i '"wl 1 fm ' » *■ ■ < 43 > 



Linear theory is expected[13] in fact to describe the evolution of the power spectrum at small k only for 
n < 4. This is because any reorganisation of smaller scales itself produces a term « k 4 , so that evolution 
on small scales is then expected to dominate the evolution of large scale modify that at small k. 



In the asymptotic behavior the relative rescaling in space for any two times becomes a 
function only of the difference in time between them so that we can write 



( r \ 2A/ 
£(r,; + A0 = £(^-^y,*J ; R s (At)=e^M. (44) 

Let us now return to Fig. 5. In addition to the measured values of R s (t) is plotted the 
best-fit (corresponding to the time t s ~ 2.5) given by Eq. (44), with n = 2 for the SL. 
The agreement is clearly very good. The system is thus indeed very well described by 
the self-similar scaling predicted by linearized fluid theory after an initial transient time. 

All these behaviours of the two point correlations are qualitatively just like those 
observed in cosmological simulations in an expanding universe. More general than 
the "self- similar" properties (which apply only to power law initial conditions), the 
clustering can be described as "hierarchical" , a feature typical of all currently favoured 
cosmological ("cold dark matter" -type) models: structures develop at a scale which 
increases in time, at a rate which can be determined from linear theory. This is given the 
following physical interpretation: clustering may be understood essentially as produced 
by the collapse of small initial overdensities which evolve as prescribed by linear theory, 
independently of pre-existing structures at smaller scales, until they "go non-linear". 
Thus the process is essentially characterised by a flow of power from large scales to 
small scales 30 . 

One of the crucial points about this description of clustering is that it tends to support 
the hypothesis that cosmological simulations may be understood almost entirely within a 
fluid description, i.e., more broadly within the framework of Vlasov-Poisson equations. 
We will discuss this point further below, after the next section, which shows that we 
can in fact understand significant aspects of the clustering we have described without 
assuming this framework. 

EVOLUTION FROM SHUFFLED LATTICE IC: THEORY 

We present in this section various analytical and numerical approaches to understanding 
more fully the dynamical evolution from SL initial conditions, beyond the linearized 
fluid theory which we have considered in the previous section. 

Phase 1: "Melting" of the lattice 

The evolution of a self-gravitating perturbed lattice — whether these perturbations 
are correlated or not — may be described using a perturbative treatment completely 
analogous to that used in condensed matter physics to describe classical phonons (see 



The extrapolation of this kind of picture has been used in constructing various quite successful phe- 
nomenological models which has been used to describe the non-linear clustering observed in numerical 
simulations, notably the so-called "halo-model" (see e.g. Ref. [26]) to which we will return below, as well 
as prescriptions like that of Peacock and Dodds [27] for calculating the final power spectrum. 



e.g. Ref. [28]). Substituting r,-(f) = R + u(R,f) directly in Eq. (1), where R is the lattice 
vector of the z'th particle and u(R, t) its displacement, one obtains, expanding each term 
in the force to linear order in the displacements, 

u(R, t) = - £ 0(R - R / )u(R / , t) . (45) 

R' 

The matrix 3> is the dynamical matrix. For gravity it may be written 

% V (R ± 0) = Gm -3^) , %v(0) = - r I %v(R) (46) 

where S^ v is the Kronecker delta. Note that a sum over the copies, associated with the 
periodic boundary conditions, is implicit in these expressions. 

The Bloch theorem for lattices tells us that @ is diagonalized by plane waves in 
reciprocal space. We define the discrete Fourier transform and its inverse by 

u(k,0 =I R e-*- R u(R,0 (47) 
u(R,0 =ilk^' k ' R fi(k,0, (48) 

where the sum in Eq. (48) is over the first Brillouin zone of the lattice, i.e., for a simple 
cubic lattice k = n(2n/L), where n is a vector of integers of which each component n ; 
(/ = 1,2,3) takes all integer values in the range —N l ^/2 < m < N 1 ^ 3 /2. Using these 
definitions in Eq. (45) we obtain 

u(M) = -#(k)u(k,f) (49) 

where ^(k), the Fourier transform (FT) of f^(R), is a symmetric 3x3 matrix for each 
k. 

Diagonalising ^(k) one can determine, for each k, three orthonormal eigenvectors 
e„(k) and their eigenvalues fi)^(k) (n = 1,2, 3), which obey the Kohn sum rule 

J>„ 2 (k) = -47rGp . (50) 

n 

Given the initial displacements and velocities at a time t = 0, the dynamical evolution 
of the particle trajectories is then given as 

U ( R >') = lX MM)fi(k,0) + J2(M)fl(k,0)] e ' kR (51) 
iV k 

where the matrix elements of the "evolution operator" & and are 

^uv(M)= lLit / n(k,0(en(k)) M (e„(k)) v (52) 
£nv(Kt)= LLiK(k,0(e„(k)) M (e„(k)) v . (53) 




and 31 

U n (k,t) = cosh ^fi)„ 2 (k)^ , ,V B (k,0 = [sinh ^0)2(k)^]/y^(k) (54) 

and the convention for the square root is: Va 2 = \ (o\ if (0 2 > 0, and \fa 2 = i\a>\ if 
(O 2 < 0. Thus, depending on the sign of C0 2 (k) and the initial conditions, each mode is 
either growing/decaying or oscillatory in nature. 

To fully solve for the evolution in this approximation requires thus only the determi- 
nation of the eigenvectors and eigenvalues of the N 3 x3 matrices J^(k). This is straight- 
forward to do and computationally inexpensive. In Fig. 7 are shown the spectrum of the 
normalized eigenvalues 

e ^-S^T' (55) 

for a 16 3 simple lattice. Full detail on these calculations may be found in Ref. [29]. 
The eigenvalues are plotted as a function of the modulus k of k, in units of the Nyquist 
frequency k^ = j. The fact that the eigenvalues at a given k do not have the same value 
is a manifestation of the anisotropy of the lattice. Shown in the plot are also lines linking 
eigenvectors oriented in some chosen directions. This allows one to see the branch 
structure of the spectrum, which is familiar in the context of analogous calculations 
in condensed matter physics (see e.g. [28]). Increasing particle number only changes 
the density of reciprocal lattice vectors in these plots (since k/k^ = n/TV 1 / 3 ), just filling 
in more densely the plot of the eigenvalues in Fig. 7 but leaving its form essentially 
unchanged. 

Formally the expansion used of the force in this treatment — which we call "particle 
linear theory" or PLT — is valid until the relative separation of any two particles 



31 When this treatment is generalized (see Ref. [29]) for cosmological simulations, only these mode 
functions change. 



becomes equal to their initial separation on the lattice, but numerical study is required 
to determine the usefulness in practice of the linearized approximation. In [29] we have 
investigated this question in detail, comparing the results of the evolution under full 
gravity (in numerical simulations), with that obtained using the formulae just given. It 
turns out that it gives an excellent description of the evolution (for the typical quantities 
describing clustering) until a time when the average relative displacement approaches 
the lattice spacing, i.e., until a significant number of particles start to approach close to 
one another (and the lattice "melts"). 



The PLT approximation for the force breaks down when particles approach closely 
nearby particles, simply because the force exerted becomes dominated by these par- 
ticles. It is thus natural to consider whether this perturbative phase may be followed 
by one in which these interactions between nearest neighbors (NN) may be relevant in 
understanding the clustering which emerges. It turns out that in fact a very good approx- 
imation to the evolution of the SL is given by assuming an abrupt switch from a phase 
in which the forces on particles are given by the PLT approximation, to a phase in each 
particle feels only the force due its NN. Further this approximation works until a very 
significant amount of non-linear correlation has developed, essentially covering most of 
the transient period prior to the asymptotic self-similar scaling we have discussed above. 

In Refs. [1,2] indirect evidence was presented for the validity of such a model, using 
the fact the following relation was observed to hold in numerical simulations, in the 
relevant range of time scales: 



Here co(r) is the NN PDF, i.e., Q)(r)dr is the probability that a particle has its NN in the 
radial shell between [r,r + dr]. The relation is exact[7] if all but two point correlations 
may be neglected, i.e., if the correlations may be fully accounted for only due to two 
body correlations. 

In a recent paper [3] we have shown directly the validity of this two phase approxima- 
tion to the evolution of an SL at early times, by comparing the full numerical evolution 
with gravity to a direct numerical integration of this model. Results are shown in Fig. 8. 
The only free parameter here is the time at which we switch from integration with PLT 
to that with the single NN, and the results in the figure are given for the choice of this 
time in each case which gives the best fit to the full evolution 32 . The time t mSLX at which 
the comparisons are given in this figure are the latest time up to which the NN approx- 
imation gives a close fit to the full evolution, i.e., the largest time up to which this two 
phase model of the evolution works well in each case. 



The values of this time determined in the way are completely consistent with the model and can be 
understood in greater detail by a closer analysis (see Ref. [3]). 



Phase 2: Nearest neighbour domination 




(56) 




FIGURE 8. The two-point correlation function at the times t max — 1,2.5,3.5,4.5 for the different SL 
initial conditions as indicated, in both the full gravity simulations (thick lines) and the simulations of the 
two phase model described in the text (thin lines). For clarity the x axis has been rescaled for each initial 
condition (as otherwise the curves are, to a very good approximation, all superimposed). From Ref. [3]. 
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FIGURE 9. Evolution of the rescaling factor R s (t) in the different simulations. Also shown is the self- 
similar behavior Eq. (43). From Ref. [1]. 



Phase 3: Self-similar evolution 

We have discussed above how the evolved SL approximates very well, at sufficiently 
long times, the "self-similar" scaling corresponding to the temporal behavior of the 
function R s (t) predicted by the linearized fluid theory. While above we showed the 
results for the SL128 simulation, we have carried out in Ref. [1] the same analysis for 
the other initial conditions, and the results are shown in Fig. 9. We see that in each case 
the tendency to join the asymptotic scaling behaviour is observed, but occurs at a later 
time. This is simply because as the initial <5 decreases (i.e. the normalised shuffling, cf. 
Table 1) the duration of the PLT phase increases, and thus the NN phase is entered later, 
and breaks down later. Comparing the estimated breakdown times t max for the two phase 
approximation given in the caption of Fig. 8 above with Fig. 9, we see that each system is 



close to joining the self-similar behavior at the time when the NN approximation breaks 
down. 

We do not have, currently, an analytical or semi-analytical treatment which allows 
us to understand the essential physical question of the origin of the functional form of 
the spatial dependence of the asymptotic correlation function. However the following 
observation, which has been emphasized in Refs. [1, 2, 3] is important. Both Figs. 9 
and 5 are derived in the approximation that the correlation function is fit by a single 
functional form. To the extent that such a fit is reasonably good, it means that the 
correlation function — in the relevant range of scale, in which £ (r) ranges from close 
to 10 2 down to about 0.1 — has approximately the spatial dependence of the asymptotic 
correlation function. This suggests strongly — but does not prove — that this asymptotic 
form of the non-linear correlations is in fact determined by the early time clustering, in 
the NN dominated phase, while its temporal behavior is determined by the fluid limit. 
We will discuss this point further in the next section, because it is of direct relevance to 
the question we consider concerning the role of discreteness in simulations of this kind. 

SOME OPEN ISSUES 

We have underlined the qualitative similarity of the evolution of the infinite SL configu- 
rations to that observed in numerical studies of gravitational clustering in the universe. 
The simulations in this context differ from those here in the initial perturbations applied 
to the lattice, and the fact that the expansion of the spatial background is incorporated 
[using the equations of motion Eq. (13) instead of Eq. (1)]. We conclude that the model 
we have studied, in which we neglect these extra complexities specific to cosmology, 
provides an interesting "toy model" to address some essential physical questions about 
this wider class of simulations. Further we have framed the model in a purely statistical 
mechanics language so that these physical issues, which are related to ones addressed in 
the wider context of long range interacting systems (as in this book), may be more easily 
accessible to non-cosmologists. 

We thus conclude with a discussion of some specific questions which we can address 
with this "toy model" and which are of direct relevance to the analogous cosmological 
simulations. We focus mostly on the issue of discreteness effects, and then briefly discuss 
the more general, but related, problem of gaining a better theoretical understanding of 
the non-linear regime of gravitational clustering in these infinite systems. 

Discreteness effects and the VP limit 

In current cosmological theories of structure formation of the universe the dominant 
clustering component is purely self-gravitating ("dark") non-relativistic ("cold") matter. 
Theoretically it is described (see e.g. [13]) by a set of Vlasov-Poisson(VP) equations. 
Although, as mentioned above, there is no rigorous derivation demonstrating the appli- 
cability of this limit in this context, it can be understood as an appropriate (mean-field) 
approximation because there is a separation of length scales, between those character- 
ising the granularity of the microscopic physical particles and the astrophysical scales 



relevant to cosmology. This allows one to view the VP equations as obtained by a coarse- 
graining over an intermediate scale (see e.g. Ref.[30]). 

Solving the VP equations for the 6-dimensional one particle phase space density is, 
however, not feasible numerically, because of the strong non-linear coupling of scales in 
the problem. For this reason cosmologists adopt the V-body approach, simulating self- 
gravitating "macro-particles" in real space. The mass of these unphysical particles are 
typically 10 60 ~ 80 larger than the supposed mass of the physical particles, and thus the 
differences between the lattice spacing i and a physical particle separation is at least of 
order 10 20 ! The problem of determining the effects of this discretization on the measured 
clustering is the "discreteness problem". Until now it is a problem which has been 
addressed mostly in a very qualitative manner, and there is considerable disagreement 
in the literature amongst the few authors who have attempted to address it (see Ref. [31] 
for references). It is now in fact a problem of very practical relevance as extremely 
precise predictions will be needed in the near future from these simulations in order 
to compare the current theoretical models with the experimental results produced by 
several different observational techniques. 

Let us just briefly consider the issue within the context of the system we have dis- 
cussed here. Is the evolution of the system described by the VP equations? While we 
have shown that the temporal dependence of the asymptotic self-similar behaviour can 
be derived within this framework, we have seen in the previous section that we can un- 
derstand aspects of the evolution using a NN approximation, which are clearly incom- 
patible with a VP limit. So what is the VP limit for our system? How do we extrapolate 
to such a limit? 

The class of infinite SL we considered is characterised, as we discussed by just two 
parameters, which we can choose to be t and 8. Gravity, however, is a scale free and 
so, dynamically, there is really only one parameter: any two SL with the same 8 are in 
fact equivalent up to a redefinition of the unit of length (which may be chosen equal 
to £). Thus 8 is in fact the single physically relevant parameter. Can we recover a VP 
approximation for some limiting value of 5? The answer is clearly in the negative: for 
any value of 8 the NN phase remains 33 . It is interesting to note also that as 8 — > the 
PLT treatment we described becomes valid for an arbitrarily long time. In fact it can be 
shown[29, 31] that this leads to a divergence of the evolution from the VP limit, given 
in this case by the linearized fluid evolution. In fact PLT approximates well the fluid 
evolution on time scales of a few dynamical times, but deviates more and more in time, 
with an asymptotic divergence 34 . 

To define a VP limit for our infinite system we clearly need a length scale. Indeed 
as we noted in our discussion of the comparison of the thermodynamic limit and the 
VP limit above, the VP limit must be one in which the particle density hq — > °o, which 
corresponds to sending £ — > 0. We must therefore have another length scale to compare 



3 We have discussed only explicitly the case 8 < 1 . For any 8 > 1 the SL at small scales is identical to a 
Poisson distribution. In this case the NN phase occurs always at the beginning of the evolution. For more 
detail, see Refs. [32, 3]. 

34 We note that this finding makes sense physically as the VP limit is taken at fixed time. We consider 
here instead the limit in which the particle density is held fixed but the time diverges. 



the density with. The scale which can evidently be used is the smoothing length £, in 
which case the VP limit is £ — > at fixed £ (and then, if the limit is indeed defined, 
£ — > 0). Studying what we have learnt about the evolution of the SL, we can verify that 
this limit does indeed lead to the recovery of the fluid limit of the two-phase model we 
have discussed. Firstly PLT, modified to include the smoothing, can be shown[29, 31] to 
converge, if £ <C £, to the fluid limit. Further the NN dominated phase we have described 
disappears also once we reach the regime £ <C £, as the force due to nearby particles is 
now no longer felt. 

The conclusion of this analysis is that to attain strictly the VP limit of such simulations 
we should work in the parameter range £ <C £. This is not the current practice in 
cosmological simulations(see e.g. Ref. [33]), but does correspond to the conclusions 
of one of the few controlled numerical studies of the issue in the literature[34]. While 
the results of (most) simulations, which use £ <C £, are not necessarily grossly wrong, 
we conclude that numerical extrapolations to the regime £ > £ will be necessary to 
determine how precise they are. 



Non-linear structures and quasi-stationary states 

We conclude with a brief discussion of the more general question of the nature of 
non-linear clustering in these systems. In recent years, in parallel with the description of 
the results of simulations in terms of the two-point correlation properties, cosmologists 
have developed a different phenomenological description of the results of simulations, 
using so-called "halo models" (see e.g. [26]). These models following naturally from the 
"hierarchical" picture of clustering, in which it is initial small overdensities at a given 
scale which evolve linearly until they collapse, essentially then behaving like indepen- 
dent finite sub-systems. The non-linear structure in simulations can be well described as 
a collection of such "halos", which are roughly spherical and approximately virialised 
structures, with properties very similar to those of structures formed after the "violent 
relaxation" of a finite gravitating system in an initially unstable configuration. These 
"halos" are believed to be understood again within the framework of VP equations, and 
are thus closely analogous to the "quasi-stationary states" (QSS) which have recently 
been much studied in certain statistical mechanical toy models (see contribution of S. 
Ruffo in this volume). While they have been characterised in great detail in huge numeri- 
cal studies, and in particular certain apparently "universal" properties identified (see e.g. 
Refs. [35, 36, 37, 38]) their physics is not at all understood. In this respect it is interest- 
ing to note that a recent detailed study [39] concludes that the framework of Lynden-Bell, 
exploited recently in the study of such QSS in simple systems, is apparently not useful 
in predicting their properties. 

In the SL system it appears visually, and would be expected from the qualitative 
features we have seen, that such a halo model description of the clustering could be 
given. It would be interesting to analyse in this context the nature of these halos, and also 
their relation to virialized structures formed in finite open systems (e.g. in spherical cold 
collapse). Further it would be perhaps useful to explore in particular lower dimensional 
(i.e. two or one) dimensional models which resemble more closely the infinite system we 



have discussed in order to simplify and possibly make closer contact with the insights 
gained, notably about QSS, in the studies of simple statistical mechanical models. 

I am indebted to my collaborators in the research reported here: Thierry Baertschiger, 
Andrea Gabrielli, Bruno Marcos and Francecso Sylos Labini. I also thank Michael 
Kiessling and Francois Sicard for several conversations which were very useful in 
relation to the discussion of the first section of this paper. 
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